home *** CD-ROM | disk | FTP | other *** search
/ Graphics Plus / Graphics Plus.iso / libs / gle / util / fitls / findmin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-11-29  |  6.2 KB  |  321 lines
  1. #include <math.h>
  2.  
  3. #define ITMAX 100
  4. #define CGOLD 0.3819660
  5. #define ZEPS 1.0e-10
  6. #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
  7. #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
  8.  
  9. double brent(ax,bx,cx,f,tol,xmin)
  10. double ax,bx,cx,tol,*xmin;
  11. double (*f)();    /* ANSI: double (*f)(double); */
  12. {
  13.     int iter;
  14.     double a,b,d,etemp,fu,fv,fw,fx,p,q,r,tol1,tol2,u,v,w,x,xm;
  15.     double e=0.0;
  16.     void nrerror();
  17.  
  18.     a=((ax < cx) ? ax : cx);
  19.     b=((ax > cx) ? ax : cx);
  20.     x=w=v=bx;
  21.     fw=fv=fx=(*f)(x);
  22.     for (iter=1;iter<=ITMAX;iter++) {
  23.         xm=0.5*(a+b);
  24.         tol2=2.0*(tol1=tol*fabs(x)+ZEPS);
  25.         if (fabs(x-xm) <= (tol2-0.5*(b-a))) {
  26.             *xmin=x;
  27.             return fx;
  28.         }
  29.         if (fabs(e) > tol1) {
  30.             r=(x-w)*(fx-fv);
  31.             q=(x-v)*(fx-fw);
  32.             p=(x-v)*q-(x-w)*r;
  33.             q=2.0*(q-r);
  34.             if (q > 0.0) p = -p;
  35.             q=fabs(q);
  36.             etemp=e;
  37.             e=d;
  38.             if (fabs(p) >= fabs(0.5*q*etemp) || 
  39.                 p <= q*(a-x) || p >= q*(b-x)) {                
  40.                 d=CGOLD*(e=(x >= xm ? a-x : b-x));
  41.             } else {
  42.                 d=p/q;
  43.                 u=x+d;
  44.                 if (u-a < tol2 || b-u < tol2)
  45.                     d=SIGN(tol1,xm-x);
  46.             }
  47.         } else {
  48.             d=CGOLD*(e=(x >= xm ? a-x : b-x));
  49.         }
  50.         u=(fabs(d) >= tol1 ? x+d : x+SIGN(tol1,d));
  51.         fu=(*f)(u);
  52.         if (fu <= fx) {
  53.             if (u >= x) a=x; else b=x;
  54.             SHFT(v,w,x,u)
  55.             SHFT(fv,fw,fx,fu)
  56.         } else {
  57.             if (u < x) a=u; else b=u;
  58.             if (fu <= fw || w == x) {
  59.                 v=w;
  60.                 w=u;
  61.                 fv=fw;
  62.                 fw=fu;
  63.             } else if (fu <= fv || v == x || v == w) {
  64.                 v=u;
  65.                 fv=fu;
  66.             }
  67.         }
  68.     }
  69.     nrerror("Too many iterations in BRENT");
  70.     *xmin=x;
  71.     return fx;
  72. }
  73.  
  74. #undef ITMAX
  75. #undef CGOLD
  76. #undef ZEPS
  77. #undef SIGN
  78. extern int ncom;    /* defined in LINMIN */
  79. extern double *pcom,*xicom,(*nrfunc)();
  80.  
  81. double f1dim(x)
  82. double x;
  83. {
  84.     int j;
  85.     double f,*xt,*vector();
  86.     void free_vector();
  87.  
  88.     xt=vector(1,ncom);
  89.     for (j=1;j<=ncom;j++) xt[j]=pcom[j]+x*xicom[j];
  90.     f=(*nrfunc)(xt);
  91.     free_vector(xt,1,ncom);
  92.     return f;
  93. }
  94. #define TOL 2.0e-4
  95.  
  96. int ncom=0;    /* defining declarations */
  97. double *pcom=0,*xicom=0,(*nrfunc)();
  98.  
  99. void linmin(p,xi,n,fret,func)
  100. double p[],xi[],*fret,(*func)();
  101. int n;
  102. {
  103.     int j;
  104.     double xx,xmin,fx,fb,fa,bx,ax;
  105.     double brent(),f1dim(),*vector();
  106.     void mnbrak(),free_vector();
  107.  
  108.     ncom=n;
  109.     pcom=vector(1,n);
  110.     xicom=vector(1,n);
  111.     nrfunc=func;
  112.     for (j=1;j<=n;j++) {
  113.         pcom[j]=p[j];
  114.         xicom[j]=xi[j];
  115.     }
  116.     ax=0.0;
  117.     xx=1.0;
  118.     bx=2.0;
  119.     mnbrak(&ax,&xx,&bx,&fa,&fx,&fb,f1dim);
  120.     *fret=brent(ax,xx,bx,f1dim,TOL,&xmin);
  121.     for (j=1;j<=n;j++) {
  122.         xi[j] *= xmin;
  123.         p[j] += xi[j];
  124.     }
  125.     free_vector(xicom,1,n);
  126.     free_vector(pcom,1,n);
  127. }
  128.  
  129. #undef TOL
  130. #include <math.h>
  131.  
  132. #define GOLD 1.618034
  133. #define GLIMIT 100.0
  134. #define TINY 1.0e-20
  135. #define MAX(a,b) ((a) > (b) ? (a) : (b))
  136. #define SIGN(a,b) ((b) > 0.0 ? fabs(a) : -fabs(a))
  137. #define SHFT(a,b,c,d) (a)=(b);(b)=(c);(c)=(d);
  138.  
  139. void mnbrak(ax,bx,cx,fa,fb,fc,func)
  140. double *ax,*bx,*cx,*fa,*fb,*fc;
  141. double (*func)();    /* ANSI: double (*func)(double); */
  142. {
  143.     double ulim,u,r,q,fu,dum;
  144.  
  145.     *fa=(*func)(*ax);
  146.     *fb=(*func)(*bx);
  147.     if (*fb > *fa) {
  148.         SHFT(dum,*ax,*bx,dum)
  149.         SHFT(dum,*fb,*fa,dum)
  150.     }
  151.     *cx=(*bx)+GOLD*(*bx-*ax);
  152.     *fc=(*func)(*cx);
  153.     while (*fb > *fc) {
  154.         r=(*bx-*ax)*(*fb-*fc);
  155.         q=(*bx-*cx)*(*fb-*fa);
  156.         u=(*bx)-((*bx-*cx)*q-(*bx-*ax)*r)/
  157.             (2.0*SIGN(MAX(fabs(q-r),TINY),q-r));
  158.         ulim=(*bx)+GLIMIT*(*cx-*bx);
  159.         if ((*bx-u)*(u-*cx) > 0.0) {
  160.             fu=(*func)(u);
  161.             if (fu < *fc) {
  162.                 *ax=(*bx);
  163.                 *bx=u;
  164.                 *fa=(*fb);
  165.                 *fb=fu;
  166.                 return;
  167.             } else if (fu > *fb) {
  168.                 *cx=u;
  169.                 *fc=fu;
  170.                 return;
  171.             }
  172.             u=(*cx)+GOLD*(*cx-*bx);
  173.             fu=(*func)(u);
  174.         } else if ((*cx-u)*(u-ulim) > 0.0) {
  175.             fu=(*func)(u);
  176.             if (fu < *fc) {
  177.                 SHFT(*bx,*cx,u,*cx+GOLD*(*cx-*bx))
  178.                 SHFT(*fb,*fc,fu,(*func)(u))
  179.             }
  180.         } else if ((u-ulim)*(ulim-*cx) >= 0.0) {
  181.             u=ulim;
  182.             fu=(*func)(u);
  183.         } else {
  184.             u=(*cx)+GOLD*(*cx-*bx);
  185.             fu=(*func)(u);
  186.         }
  187.         SHFT(*ax,*bx,*cx,u)
  188.         SHFT(*fa,*fb,*fc,fu)
  189.     }
  190. }
  191.  
  192. #undef GOLD
  193. #undef GLIMIT
  194. #undef TINY
  195. #undef MAX
  196. #undef SIGN
  197. #undef SHFT
  198. #include <math.h>
  199.  
  200. #define ITMAX 200
  201. static double sqrarg;
  202. #define SQR(a) (sqrarg=(a),sqrarg*sqrarg)
  203.  
  204. void powell(double p[],double **xi, int n, double ftol, int *iter, double *fret, double (*func)());
  205. void powell(p,xi,n,ftol,iter,fret,func)
  206. double p[],**xi,ftol,*fret,(*func)();
  207. int n,*iter;
  208. {
  209.     int i,ibig,j;
  210.     double t,fptt,fp,del;
  211.     double *pt,*ptt,*xit,*vector();
  212.     void linmin(),nrerror(),free_vector();
  213.  
  214.     pt=vector(1,n);
  215.     ptt=vector(1,n);
  216.     xit=vector(1,n);
  217.     *fret=(*func)(p);
  218.     for (j=1;j<=n;j++) pt[j]=p[j];
  219.     for (*iter=1;;(*iter)++) {
  220.         fp=(*fret);
  221.         ibig=0;
  222.         del=0.0;
  223.         for (i=1;i<=n;i++) {
  224.             for (j=1;j<=n;j++) xit[j]=xi[j][i];
  225.             fptt=(*fret);
  226.             linmin(p,xit,n,fret,func);
  227.             if (fabs(fptt-(*fret)) > del) {
  228.                 del=fabs(fptt-(*fret));
  229.                 ibig=i;
  230.             }
  231.         }
  232.         if (2.0*fabs(fp-(*fret)) <= ftol*(fabs(fp)+fabs(*fret))) {
  233.             free_vector(xit,1,n);
  234.             free_vector(ptt,1,n);
  235.             free_vector(pt,1,n);
  236.             return;
  237.         }
  238.         if (*iter == ITMAX) nrerror("Too many iterations in routine POWELL");
  239.         for (j=1;j<=n;j++) {
  240.             ptt[j]=2.0*p[j]-pt[j];
  241.             xit[j]=p[j]-pt[j];
  242.             pt[j]=p[j];
  243.         }
  244.         fptt=(*func)(ptt);
  245.         if (fptt < fp) {
  246.             t=2.0*(fp-2.0*(*fret)+fptt)*SQR(fp-(*fret)-del)-del*SQR(fp-fptt);
  247.             if (t < 0.0) {
  248.                 linmin(p,xit,n,fret,func);
  249.                 for (j=1;j<=n;j++) xi[j][ibig]=xit[j];
  250.             }
  251.         }
  252.     }
  253. }
  254.  
  255. #undef ITMAX
  256. #undef SQR
  257.  
  258. #ifdef __TURBOC__
  259. #include <alloc.h>
  260. #endif
  261. #include <stdio.h>
  262.  
  263. void nrerror(error_text)
  264. char error_text[];
  265. {
  266.     void exit();
  267.  
  268.     fprintf(stderr,"Numerical Recipes run-time error...\n");
  269.     fprintf(stderr,"%s\n",error_text);
  270.     fprintf(stderr,"...now exiting to system...\n");
  271.     exit(1);
  272. }
  273.  
  274.  
  275.  
  276. double *vector(nl,nh)
  277. int nl,nh;
  278. {
  279.     double *v;
  280.  
  281.     v=(double *)malloc((unsigned) (nh-nl+1)*sizeof(double));
  282.     if (!v) nrerror("allocation failure in vector()");
  283.     return v-nl;
  284. }
  285.  
  286. double **matrix(nrl,nrh,ncl,nch)
  287. int nrl,nrh,ncl,nch;
  288. {
  289.     int i;
  290.     double **m;
  291.  
  292.     m=(double **) malloc((unsigned) (nrh-nrl+1)*sizeof(double*));
  293.     if (!m) nrerror("allocation failure 1 in matrix()");
  294.     m -= nrl;
  295.  
  296.     for(i=nrl;i<=nrh;i++) {
  297.         m[i]=(double *) malloc((unsigned) (nch-ncl+1)*sizeof(double));
  298.         if (!m[i]) nrerror("allocation failure 2 in matrix()");
  299.         m[i] -= ncl;
  300.     }
  301.     return m;
  302. }
  303.  
  304. void free_vector(v,nl,nh)
  305. double *v;
  306. int nl,nh;
  307. {
  308.     free((char*) (v+nl));
  309. }
  310.  
  311. void free_matrix(m,nrl,nrh,ncl,nch)
  312. double **m;
  313. int nrl,nrh,ncl,nch;
  314. {
  315.     int i;
  316.  
  317.     for(i=nrh;i>=nrl;i--) free((char*) (m[i]+ncl));
  318.     free((char*) (m+nrl));
  319. }
  320.                                                                                                                                                        
  321.